####### Main results ##############

# This script produces:
# Table 1 (main regression results)
# Figure 4
# Figure 5
# Table 2 (Difference-in-differences results)



# Session info
# R version 4.1.0 (2021-05-18)
# Platform: x86_64-apple-darwin17.0 (64-bit)
# Running under: macOS 13.1


# Load data ---------------------------------------------------------------

# Define protest wave threshold
threshold = 8
load(paste0("data/navco_week_replication",threshold,".Rda"))


# add variables for Difference-in-differences analysis
nav_did <- navco_week %>% 
  group_by(move_id) %>% 
  mutate(con_first = ifelse(row_number() == 1, con_dummy, cumsum(con_dummy) + lag(con_dummy)),
         con_first = ifelse(con_first == 1, 1, 0),
         t_week = ifelse(sum(con_dummy) > 0, week[con_first == 1], 10000),
         treat = ifelse(sum(con_first) > 0, 1, 0)) %>% 
  ungroup() %>% 
  mutate(time_to_treat =  week - t_week,
         time_to_treat = ifelse(time_to_treat < -1000, 0, time_to_treat),
         post = ifelse(time_to_treat > -1, 1, 0),
         post_t = post*treat)



# Table 1 ---------------------
# Negative binomial models 
# At the grid level
# Protest wave threshold: 7 days between protests


# run negative binomial models
  nb1 <- glm.nb(protests_future ~ con_sum,
                data=navco_week)
  
  nb2 <- glm.nb(protests_future ~ con_sum +
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence + 
                  gid_f + year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  
  nb3 <- glm.nb(protests_future ~ con_sum + 
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence + 
                  #time
                  diff +diff2  + protests_past + pr_nr_c + 
                  v2x_clphy + v2x_polyarchy +couprisk + nlights_calib_mean + gdp_log + pop_log +
                  gid_f + year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  nb4 <- glm.nb(protests_future ~ con_sum +
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence + 
                  #time
                  diff + diff2  +protests_past + pr_nr_c + 
                  factor(gid_year),
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))  
  
  
# Cluster standard errors by grids
  clustered_se <- list(
    sqrt(diag(vcovCL(nb1, vcov = vcovCL(nb1, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb2, vcov = vcovCL(nb2, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb3, vcov = vcovCL(nb3, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb4, vcov = vcovCL(nb4, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T)))
  )
  
 # Regression table
  stargazer(nb1, nb2,nb3,nb4,
            se = clustered_se,
            type = "latex",
            float = F,
            float.env = "table",
            no.space = T,
            dep.var.caption =paste0("Sum of protests in the next week"),
            dep.var.labels = c(""),
            model.numbers = F,
            column.labels = c("1", "2", "3", "4"),
            title = paste0("Concessions and protests in the next 7 days. Negative binomial models"), # is not written in latex if float=F
            notes.append = T,
            notes = c("Standard errors are clustered at the level of grid cells."),
            covariate.labels = c("Concession (sum)",
                                 "Number of protesters (sum, log)",
                                 "Protester violence (sum)",
                                 "Demand: Policy change",
                                 "Demand: Institutional reforms",
                                 "Demand: Regime change",
                                 "State violence (sum)",
                                 "Days since last protest",
                                 "Days since last protest (squared)",
                                 "Protest last week",
                                 "Protests in the rest of the country",
                                 "Physical violence index",
                                 "Democracy index",
                                 "Coup threat",
                                 "Nightlights",
                                 "GDP pc (lag and log)",
                                 "Population (lag and log)"),
            omit= c('year', 'gid', 'country'),
            omit.stat = c("logrank", "wald", "ll", "theta"),
            add.lines = list(c("Grid dummies", "No", "Yes","Yes", "Yes"),
                             c("Year dummies", "No", "Yes","Yes", "Yes"),
                             c("Grid-year interaction", "No", "No", "No", "Yes"),
                             c("Number of grids", "169", "169", "169", "169")),
            label = "table:main-negbin-national",
            out = "output/tables/Table1.tex")
  

# Figure 4 --------------------------------------------------

  # re-run model 2 
  nb2 <- glm.nb(protests_future ~ con_sum +
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence + 
                  gid_f + year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  # marginal effects plot
  plot_cap(nb2, condition = c("con_sum")) +
    labs(x = "Concessions per week", y = "Predicted # of protests in the next 7 days") +
    scale_x_continuous(breaks = seq(0,8,1)) +
    scale_y_continuous(breaks = seq(0,12,2)) +
    theme_lh
  
  ggsave("output/figures/Figure4.pdf", scale = 1,
         dpi = "retina",
         width = 7,
         height = 5)
  

# Table 2 -----------------------------------------------------------------
# Difference-in-differences models
  
  mdid1 = feols(protests_future ~ post_t |  
                  gid  + week_calendar,                             
                vcov = ~gid,                          
                data = nav_did)

  mdid2 = feols(protests_future ~ post_t +
                  # protest characteristic
                  size_log + protesterviolence + goal_policy + goal_reform + goal_regime + state_violence +
                  # previous dynamics
                  pr_nr_c + diff + diff2 + protests_past|  
                  gid + week_calendar,                             
                vcov = ~gid,  
                data = nav_did)
  
 
  etable(mdid1, mdid2,  tex = T, style.tex = style.tex("aer"),
         title = "Concessions and sum of protests in the next week (7 days). Difference-in-differences specification",
         drop = c("size_log", "protesterviolence", "goal_policy", "goal_reform", "goal_regime",
                  "state_violence", "pr_nr_c", "diff", "diff2", "protests_past", "gid", "year"),
         digits = 2,
         notes = "Standard errors are clustered at the level of grid cells.",
         fixef_sizes = T,
         fitstat = ~ r2 + n)
  
  
  #Further processing of the style of the table was done in Latex
  # Here the final latex code used in the analysis
  
  # \begin{table}[htbp]
  # \caption{Concessions and sum of protests in the next week (7 days). \\ Difference-in-differences specification}
  # \label{tab:did}
  # \bigskip
  # \centering
  # \begin{tabular}{lcc}
  # \toprule
  # & \multicolumn{2}{c}{Sum of protests in the next week}\\
  # & (1)  & (2)\\ 
  # \midrule 
  # Concession $\times$ Post & 1.04$^{**}$ & 0.38$^{*}$\\ 
  # & (0.48) & (0.21)\\ 
  # \\ \midrule
  # Grid FE (169) & Yes & Yes\\ 
  # Week FE (525) & Yes & Yes\\
  # Protest characteristics & No & Yes \\
  # Previous dynamics & No & Yes \\
  # R$^2$  & 0.66 & 0.74\\ 
  # Observations & 1,495  & 1,495\\ 
  # \bottomrule
  # \end{tabular}
  # \par \raggedright 
  # \footnotesize \textit{Note:} Standard errors are clustered at the level of grid cells. \\
  # {$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01} 
  # \end{table}
  # 
  
  
  
